
set more off 
cap log close 
clear all
set niceness 10 
set processor 4

*** Paths ****
gl path1 = "/Users/Wei/Dropbox/Census"
gl path2 = "/Users/Wei/Dropbox/Fertility/Workingdata"
gl path3 = "/Users/Wei/Dropbox/Fertility/Results"
gl path4 = "/Users/Wei/Dropbox/Fertility/Figures"

log using "$path3/Robust.log", replace 

* Figure c4 

**************************
use  "$path2/marr_policy", clear
gen early_marr = 1 - late_marr
replace high_occ = 0 if work == 0
su senior college married_ever early_marr work high_occ [aw = wt] if men, sep(100)
su senior college married_ever early_marr work high_occ [aw = wt] if women, sep(100)
su senior college married_ever early_marr work high_occ [aw = wt] if men & age <= 30, sep(100)
su senior college married_ever early_marr work high_occ [aw = wt] if men & age > 30, sep(100)
su senior college married_ever early_marr work high_occ [aw = wt] if women & age <= 30, sep(100)
su senior college married_ever early_marr work high_occ [aw = wt] if women & age > 30, sep(100)

use "$path2/fines_economic", clear 
drop lngdp_pc lngdp secondary primary
merge 1:1 prov year using "$path2/60years", keepusing(gdp_pc teacher_primary teacher_second)
ren teacher_primary primary 
ren teacher_second secondary
gen lngdp_pc = ln(gdp_pc)
save "$path2/fines_economic_new",replace 

use "$path2/marr_policy", clear
collapse men, by(prov year_birth)
gen sex_ratio = men/(1-men)
ren men male_prop
save "$path2/sex_ratio_prov_year", replace 



set more off

use "$path2/marr_policy", clear
gen mort_rate = 100 -sur_rate if women == 1 
ren year yr_t
gen year = year_birth +15
merge m:1 prov year using "$path2/fines_economic_new", keepusing(lngdp_pc secondary primary)
drop if _merge == 2 
drop _merge 
drop year
ren yr_t year
egen male_prop = mean(men), by(prov year_birth)
gen sex_ratio = male_prop/(1-male_prop)


egen fine_6_20 = rowmean(fine_age6-fine_age20)
drop fine_age*
gl FINE_ALL = "fine_6_15 fine_16_20 fine_21_25 fine_26_30"
gl FINE_VAR_4 =  "fine_6_15 fine_16_20 fine_21_30"

gl FINE_VAR_3 =  "fine_6_15 fine_16_20 fine_21_25"
gl FINE_VAR_1 = "fine_6_20"

replace nbroth = 99 if mi(nbroth)
replace nsister = 99 if mi(nsister)

gen lnsec = ln(secondary)
gen lnpri = ln(primary)
gen lnsec_mi = mi(lnsec)
gen lnpri_mi = mi(lnpri)
gen lngdp_pc_mi = mi(lnpri)

egen mean_lnsec = mean(lnsec)
egen mean_lnpri = mean(lnpri)
egen mean_lngdp = mean(lngdp_pc)

replace lnsec = mean_lnsec if mi(lnsec)
replace lnpri = mean_lnpri if mi(lnpri)
replace lngdp_pc = mean_lngdp if mi(lngdp_pc)
drop mean_lnsec  mean_lnpri mean_lngdp

gen provid = prov 
merge m:1 provid using "$path2/Comp_Law", keepusing(law_*) nogen 

cap drop elig
gen elig = 1 if year_birth >= law_year - 6 
replace elig = 0 if year_birth <= law_year - 16
replace elig = (16-(law_year - year_birth ))/10 if law_year - 16<year_birth & year_birth<law_year - 6 
drop provid


gl CONTROL_1 = "women#prov##c.year_birth women#year#year_birth"
gl CONTROL_2 = "women#prov##c.year_birth women#year#year_birth women#prov#year"

gen urban_hk = hktype == 2 
gen urban_res = city_flag <= 2 if !mi(city_flag)
replace urban_res = urban_hk if  mi(urban_res)
set more off 
replace high_occ = 0 if high_occ == . & !mi(work)

egen han_p = mean(han), by(year hhid prov region)
replace late_marr = 1- late_marr // early married 

set more off
gen nsib = nbroth + nsister
cap erase "$path3/OCP_women_census_control.xls"
cap erase "$path3/OCP_women_census_control.txt"
foreach y_var in "senior" "college ""married_ever""late_marr" "work" "high_occ"  {
local j = 1 
gen coef_men_`y_var' = .
gen coef_women_`y_var' = .
gen se_men_`y_var' = .
gen se_women_`y_var' = .
gen mean_`y_var'  = . 
gen std_`y_var'  = . 
su `y_var'  [aw = wt] if han 
replace mean_`y_var' = `r(mean)'  
replace std_`y_var' = `r(sd)'  
local j = 1 
foreach x_var in  "$FINE_VAR_1"{
reghdfe `y_var' women#c.(`x_var') if han [aw = wt], a($CONTROL_1) cluster(prov)
replace coef_men_`y_var' = _b[0.women#`x_var']  in `j'
replace coef_women_`y_var' = _b[1.women#`x_var']   in `j'
replace se_men_`y_var' = _se[0.women#`x_var']  in `j'
replace se_women_`y_var' = _se[1.women#`x_var']   in `j'
outreg2 using "$path3/OCP_women_census_control.xls", dec(3) append
local ++j 

reghdfe `y_var' women#c.(`x_var') if han  [aw = wt], a($CONTROL_1 women#prov#year) cluster(prov)
replace coef_men_`y_var' = _b[0.women#`x_var']  in `j'
replace coef_women_`y_var' = _b[1.women#`x_var'] in `j'
replace se_men_`y_var' = _se[0.women#`x_var']  in `j'
replace se_women_`y_var' = _se[1.women#`x_var']   in `j'
outreg2 using "$path3/OCP_women_census_control.xls", dec(3) append
local ++j 

reghdfe `y_var' women#c.(`x_var') if han  [aw = wt], a(women#nsib $CONTROL_1) cluster(prov)
replace coef_men_`y_var' = _b[0.women#`x_var'] in `j'
replace coef_women_`y_var' = _b[1.women#`x_var']  in `j'
replace se_men_`y_var' = _se[0.women#`x_var'] in `j'
replace se_women_`y_var' = _se[1.women#`x_var']  in `j'
outreg2 using "$path3/OCP_women_census_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' sex_ratio) if han  [aw = wt], a($CONTROL_1) cluster(prov)
replace coef_men_`y_var' = _b[0.women#`x_var'] in `j'
replace coef_women_`y_var' = _b[1.women#`x_var']  in `j'
replace se_men_`y_var' = _se[0.women#`x_var']  in `j'
replace se_women_`y_var' = _se[1.women#`x_var']  in `j'
outreg2 using "$path3/OCP_women_census_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' lngdp_pc) if han  [aw = wt], a(women#lngdp_pc_mi $CONTROL_1) cluster(prov)
replace coef_men_`y_var' = _b[0.women#`x_var']  in `j'
replace coef_women_`y_var' = _b[1.women#`x_var']   in `j'
replace se_men_`y_var' = _se[0.women#`x_var']  in `j'
replace se_women_`y_var' = _se[1.women#`x_var']   in `j'
outreg2 using "$path3/OCP_women_census_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' lnsec lnpri) if han  [aw = wt], a(women#lnsec_mi  women#lnpri_mi women#nsib $CONTROL_1) cluster(prov)
replace coef_men_`y_var' = _b[0.women#`x_var']  in `j'
replace coef_women_`y_var' = _b[1.women#`x_var']  in `j'
replace se_men_`y_var' = _se[0.women#`x_var']  in `j'
replace se_women_`y_var' = _se[1.women#`x_var']  in `j'
outreg2 using "$path3/OCP_women_census_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' elig) if han  [aw = wt], a( $CONTROL_1) cluster(prov)
replace coef_men_`y_var' = _b[0.women#`x_var'] in `j'
replace coef_women_`y_var' = _b[1.women#`x_var']  in `j'
replace se_men_`y_var' = _se[0.women#`x_var']  in `j'
replace se_women_`y_var' = _se[1.women#`x_var']  in `j'
outreg2 using "$path3/OCP_women_census_control.xls", dec(3) append
local ++j


reghdfe `y_var' women#c.(`x_var' sex_ratio lngdp_pc lnsec lnpri elig) if han [aw = wt], a(women#lngdp_pc_mi women#lnsec_mi  women#lnpri_mi  women#nsib $CONTROL_1) cluster(prov)
replace coef_men_`y_var' = _b[0.women#`x_var']  in `j'
replace coef_women_`y_var' = _b[1.women#`x_var']   in `j'
replace se_men_`y_var' = _se[0.women#`x_var']  in `j'
replace se_women_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_women_census_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' sex_ratio lngdp_pc lnsec lnpri elig) if han [aw = wt], a(women#lngdp_pc_mi women#lnsec_mi  women#lnpri_mi  women#nsib $CONTROL_1 women#prov#year) cluster(prov)
replace coef_men_`y_var' = _b[0.women#`x_var']   in `j'
replace coef_women_`y_var' = _b[1.women#`x_var']   in `j'
replace se_men_`y_var' = _se[0.women#`x_var']  in `j'
replace se_women_`y_var' = _se[1.women#`x_var']   in `j'
outreg2 using "$path3/OCP_women_census_control.xls", dec(3) append
local ++j

}
}

keep coef_* se_* mean_* std_*  
drop if mi(coef_men_senior)
gen pos = _n

save "$path3/coef_census_ind_control",replace 
 

set more off

use "$path2/marr_policy", clear
gen mort_rate = 100 -sur_rate if women == 1 
ren year yr_t
gen year = year_birth +15
merge m:1 prov year using "$path2/fines_economic_new", keepusing(lngdp_pc secondary primary)
drop if _merge == 2 
drop _merge 
drop year
ren yr_t year
egen male_prop = mean(men), by(prov year_birth)
gen sex_ratio = male_prop/(1-male_prop)


egen fine_6_20 = rowmean(fine_age6-fine_age20)
drop fine_age*
gl FINE_ALL = "fine_6_15 fine_16_20 fine_21_25 fine_26_30"
gl FINE_VAR_4 =  "fine_6_15 fine_16_20 fine_21_30"

gl FINE_VAR_3 =  "fine_6_15 fine_16_20 fine_21_25"
gl FINE_VAR_1 = "fine_6_20"

replace nbroth = 99 if mi(nbroth)
replace nsister = 99 if mi(nsister)

gen lnsec = ln(secondary)
gen lnpri = ln(primary)
gen lnsec_mi = mi(lnsec)
gen lnpri_mi = mi(lnpri)
gen lngdp_pc_mi = mi(lnpri)

egen mean_lnsec = mean(lnsec)
egen mean_lnpri = mean(lnpri)
egen mean_lngdp = mean(lngdp_pc)

replace lnsec = mean_lnsec if mi(lnsec)
replace lnpri = mean_lnpri if mi(lnpri)
replace lngdp_pc = mean_lngdp if mi(lngdp_pc)
drop mean_lnsec  mean_lnpri mean_lngdp

gen provid = prov 
merge m:1 provid using "$path2/Comp_Law", keepusing(law_*) nogen 

cap drop elig
gen elig = 1 if year_birth >= law_year - 6 
replace elig = 0 if year_birth <= law_year - 16
replace elig = (16-(law_year - year_birth ))/10 if law_year - 16<year_birth & year_birth<law_year - 6 
drop provid


gl CONTROL_1 = "women#prov##c.year_birth women#year#year_birth"
gl CONTROL_2 = "women#prov##c.year_birth women#year#year_birth women#prov#year"

gen urban_hk = hktype == 2 
gen urban_res = city_flag <= 2 if !mi(city_flag)
replace urban_res = urban_hk if  mi(urban_res)
set more off 
replace high_occ = 0 if high_occ == . & !mi(work)

egen han_p = mean(han), by(year hhid prov region)


drop if han_p >0 & han_p < 1 // drop mixed households
drop if women == 0 
set more off
 
set more off
gen nsib = nbroth + nsister
gen n_death = n_birth - n_survive

replace n_death = . if n_death < 0
replace n_birth = 0 if married_ever == 0 
replace n_mbirth = 0 if married_ever == 0  
replace n_fbirth = 0 if married_ever == 0  
replace n_survive = 0 if married_ever == 0 
replace n_mchild = 0 if married_ever == 0 
replace n_fchild = 0 if married_ever == 0 
drop if mi(n_birth)
drop if mi(n_death)


cap erase "$path3/OCP_women_birth_control.xls"
cap erase "$path3/OCP_women_birth_control.txt"
foreach y_var in "n_birth" {
gen coef_women_`y_var' = .
gen se_women_`y_var' = .
gen mean_`y_var'  = . 
gen std_`y_var'  = . 
su `y_var'  [aw = wt] if han 
replace mean_`y_var' = `r(mean)'  
replace std_`y_var' = `r(sd)'  
local j = 1 

foreach x_var in  "$FINE_VAR_1"{
reghdfe `y_var' women#c.(`x_var') if han [aw = wt], a($CONTROL_1) cluster(prov)
replace coef_women_`y_var' = _b[1.women#`x_var']   in `j'
replace se_women_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_women_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var') if han  [aw = wt], a($CONTROL_1 women#prov#year) cluster(prov)
replace coef_women_`y_var' = _b[1.women#`x_var']    in `j'
replace se_women_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_women_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var') if han  [aw = wt], a(women#nsib $CONTROL_1) cluster(prov)
replace coef_women_`y_var' = _b[1.women#`x_var']   in `j'
replace se_women_`y_var' = _se[1.women#`x_var']   in `j'
outreg2 using "$path3/OCP_women_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' sex_ratio) if han  [aw = wt], a($CONTROL_1) cluster(prov)
replace coef_women_`y_var' = _b[1.women#`x_var']   in `j'
replace se_women_`y_var' = _se[1.women#`x_var']   in `j'
outreg2 using "$path3/OCP_women_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' lngdp_pc) if han  [aw = wt], a(women#lngdp_pc_mi $CONTROL_1) cluster(prov)
replace coef_women_`y_var' = _b[1.women#`x_var']    in `j'
replace se_women_`y_var' = _se[1.women#`x_var']  in `j'
outreg2 using "$path3/OCP_women_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' lnsec lnpri) if han  [aw = wt], a(women#lnsec_mi  women#lnpri_mi  $CONTROL_1) cluster(prov)
replace coef_women_`y_var' = _b[1.women#`x_var']   in `j'
replace se_women_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_women_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' elig) if han  [aw = wt], a( $CONTROL_1) cluster(prov)
replace coef_women_`y_var' = _b[1.women#`x_var']    in `j'
replace se_women_`y_var' = _se[1.women#`x_var']   in `j'
outreg2 using "$path3/OCP_women_birth_control.xls", dec(3) append
local ++j


reghdfe `y_var' women#c.(`x_var' lngdp_pc lnsec lnpri elig sex_ratio) if han [aw = wt], a(women#lnsec_mi  women#lnpri_mi women#lngdp_pc_mi women#nsib $CONTROL_1 women#prov#year) cluster(prov)
replace coef_women_`y_var' = _b[1.women#`x_var']   in `j'
replace se_women_`y_var' = _se[1.women#`x_var']  in `j'
outreg2 using "$path3/OCP_women_birth_control.xls", dec(3) append
local ++j
}
}


keep coef_* se_* mean_* std_*  
drop if mi(coef_women_n_birth)
gen pos = _n

save "$path3/coef_census_birth_control",replace 
 
*** Children outcomes ****
use "$path2/child_gender_mort", clear 
keep if han_p == 1 
gl CONTROL_1 = "prov##c.year_birth year#year_birth"
gl CONTROL_2 = "prov##c.year_birth year#year_birth prov#year"
gen wt = 1 
replace wt = 4 if year == 2005
replace die = die*100

ren year yr_t
gen year = year_birth +15
merge m:1 prov year using "$path2/fines_economic_new", keepusing(lngdp_pc secondary primary)
drop if _merge == 2 
drop _merge 
drop year
ren yr_t year
merge m:1 prov year_birth using "$path2/sex_ratio_prov_year", keepusing(male_prop sex_ratio) 
drop if _merge == 2 
drop _merge 

replace nbroth = 99 if mi(nbroth)
replace nsister = 99 if mi(nsister)

gen lnsec = ln(secondary)
gen lnpri = ln(primary)
gen lnsec_mi = mi(lnsec)
gen lnpri_mi = mi(lnpri)
gen lngdp_pc_mi = mi(lnpri)

egen mean_lnsec = mean(lnsec)
egen mean_lnpri = mean(lnpri)
egen mean_lngdp = mean(lngdp_pc)

replace lnsec = mean_lnsec if mi(lnsec)
replace lnpri = mean_lnpri if mi(lnpri)
replace lngdp_pc = mean_lngdp if mi(lngdp_pc)
drop mean_lnsec  mean_lnpri mean_lngdp

gen provid = prov 
merge m:1 provid using "$path2/Comp_Law", keepusing(law_*) nogen 

cap drop elig
gen elig = 1 if year_birth >= law_year - 6 
replace elig = 0 if year_birth <= law_year - 16
replace elig = (16-(law_year - year_birth ))/10 if law_year - 16<year_birth & year_birth<law_year - 6 
drop provid
gen women = 1
gen nsib = nbroth + nsister

cap erase "$path3/OCP_child_birth_control.xls"
cap erase "$path3/OCP_child_birth_control.txt"
foreach y_var in "die" {
cap gen coef_child_`y_var' = .
cap gen se_child_`y_var' = .
gen mean_`y_var'  = . 
gen std_`y_var'  = . 
su `y_var'  [aw = wt] if han 
replace mean_`y_var' = `r(mean)'  
replace std_`y_var' = `r(sd)'  
local j = 1 

foreach x_var in  "$FINE_VAR_1"{
reghdfe `y_var' women#c.(`x_var') if han [aw = wt], a($CONTROL_1) cluster(prov)
replace coef_child_`y_var' = _b[1.women#`x_var']   in `j'
replace se_child_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_child_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var') if han  [aw = wt], a($CONTROL_1 women#prov#year) cluster(prov)
replace coef_child_`y_var' = _b[1.women#`x_var']   in `j'
replace se_child_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_child_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var') if han  [aw = wt], a(women#nsib $CONTROL_1) cluster(prov)
replace coef_child_`y_var' = _b[1.women#`x_var']   in `j'
replace se_child_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_child_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' sex_ratio) if han  [aw = wt], a($CONTROL_1) cluster(prov)
replace coef_child_`y_var' = _b[1.women#`x_var']   in `j'
replace se_child_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_child_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' lngdp_pc) if han  [aw = wt], a(women#lngdp_pc_mi $CONTROL_1) cluster(prov)
replace coef_child_`y_var' = _b[1.women#`x_var']   in `j'
replace se_child_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_child_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' lnsec lnpri) if han  [aw = wt], a(women#lnsec_mi  women#lnpri_mi  $CONTROL_1) cluster(prov)
replace coef_child_`y_var' = _b[1.women#`x_var']   in `j'
replace se_child_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_child_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' elig) if han  [aw = wt], a(  $CONTROL_1) cluster(prov)
replace coef_child_`y_var' = _b[1.women#`x_var']   in `j'
replace se_child_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_child_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' lngdp_pc lnsec lnpri elig sex_ratio) if han [aw = wt], a(women#lnsec_mi  women#lnpri_mi women#lngdp_pc_mi women#nsib $CONTROL_1) cluster(prov)
replace coef_child_`y_var' = _b[1.women#`x_var']   in `j'
replace se_child_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_child_birth_control.xls", dec(3) append
local ++j

reghdfe `y_var' women#c.(`x_var' lngdp_pc lnsec lnpri elig sex_ratio) if han [aw = wt], a(women#lnsec_mi  women#lnpri_mi women#lngdp_pc_mi women#nsib $CONTROL_1 women#prov#year) cluster(prov)
replace coef_child_`y_var' = _b[1.women#`x_var']   in `j'
replace se_child_`y_var' = _se[1.women#`x_var']    in `j'
outreg2 using "$path3/OCP_child_birth_control.xls", dec(3) append
local ++j
}
}
egen n_birth = sum(1), by(hhid)

keep coef_* se_* mean_* std_*  
drop if mi(coef_child_die)
gen pos = _n

save "$path3/coef_child_birth_control",replace 


* Robustness to adding more controls
use "$path2/marr_policy", clear
foreach y_var in "senior" "college" "married_ever" "late_marr" "work" "high_occ"{
su `y_var' [aw = wt] if men & han
sca std_men_`y_var' = `r(sd)'
su `y_var' [aw = wt] if women & han
sca std_women_`y_var' = `r(sd)'
}
sca list
use "$path3/coef_census_ind_control", clear
foreach y_var in "senior" "college" "married_ever" "late_marr" "work" "high_occ"{
replace coef_men_`y_var' = coef_men_`y_var'/std_men_`y_var'
replace se_men_`y_var' = se_men_`y_var'/std_men_`y_var'
replace coef_women_`y_var' = coef_women_`y_var'/std_women_`y_var'
replace se_women_`y_var' = se_women_`y_var'/std_women_`y_var'

gen up_men_`y_var' = coef_men_`y_var' + 1.65*se_men_`y_var'
gen low_men_`y_var' = coef_men_`y_var' - 1.65*se_men_`y_var'
gen up_women_`y_var' = coef_women_`y_var' + 1.65*se_women_`y_var'
gen low_women_`y_var' = coef_women_`y_var' - 1.65*se_women_`y_var'
}
drop if p == 8 
replace p = _n
gen label = ""
replace label = "Original" if pos == 1 
replace label = "Prov-Year FE controlled" if pos == 2
replace label = "# of siblings controlled" if pos == 3 
replace label = "Sex ratio at birth controlled" if pos == 4
replace label = "GDP /capita at age 15 controlled" if pos == 5 
replace label = "Education resources controlled" if pos == 6 
replace label = "CSLs controlled"  if pos == 7
replace label = "All above" if pos == 8 


replace p = -p
foreach gender in "women"{
tw (scatter p coef_`gender'_senior , xscale(alt) xline(0) xtit("Coef. on Fines at 6-20 / Std", size(small))  mlabel(label) mlabp(12)  xlabel(-0.2(0.1)0.4, grid  )) ///
(rcap up_`gender'_senior low_`gender'_senior p  , ytit(" ") horizontal lp(dash) lc(green)  legend(pos(5) ring(0) col(1) size(small) lab(1 "Senior high school") lab(2 "90% CI"))  ylabel(-9 " " 0 " ") )
gr export "$path4/fig_c4a.eps",replace

tw (scatter p coef_`gender'_college , xscale(alt) xline(0) xtit("Coef. on Fines at 6-20 / Std", size(small))  mlabel(label) mlabp(12)  xlabel(-0.2(0.1)0.4, grid)) ///
(rcap up_`gender'_college low_`gender'_college p  , ytit(" ") horizontal lp(dash) lc(green)  legend(pos(5) ring(0) col(1) size(small) lab(1 "College") lab(2 "90% CI"))  ylabel(-9 " " 0 " ") )
gr export "$path4/fig_c4b.eps",replace


tw (scatter p coef_`gender'_late_marr, xscale(alt) xline(0) xtit("Coef. on Fines at 6-20 / Std", size(small))  mlabel(label) mlabp(12)  xlabel(-0.4(0.2)1, grid )) ///
(rcap up_`gender'_late_marr low_`gender'_late_marr p  , ytit(" ") horizontal lp(dash) lc(green)  legend(pos(5) ring(0) col(1) size(small) lab(1 "Married before 25") lab(2 "90% CI"))  ylabel(-9 " " 0 " ") )
gr export "$path4/fig_c4c.eps",replace


tw (scatter p coef_`gender'_high_occ, xscale(alt) xline(0) xtit("Coef. on Fines at 6-20 / Std", size(small))  mlabel(label) mlabp(12)  xlabel(-0.2(0.1)0.4, grid )) ///
(rcap up_`gender'_high_occ low_`gender'_high_occ p  , ytit(" ") horizontal lp(dash) lc(green)  legend(pos(5) ring(0) col(1) size(small) lab(1 "White-collar") lab(2 "90% CI"))  ylabel(-9 " " 0 " ") )
gr export "$path4/fig_c4d.eps",replace
}



use "$path2/marr_policy", clear
keep if women == 1 & married_ever == 1
gen mort_rate = 100 - sur_rate
foreach y_var in "n_birth" {
su `y_var' [aw = wt] if women & han
sca std_women_`y_var' = `r(sd)'
}
sca list

use "$path3/coef_census_birth_control", clear
foreach y_var in "n_birth" {
replace coef_women_`y_var' = coef_women_`y_var'/std_`y_var'
replace se_women_`y_var' = se_women_`y_var'/std_`y_var'
gen up_women_`y_var' = coef_women_`y_var' + 1.65*se_women_`y_var'
gen low_women_`y_var' = coef_women_`y_var' - 1.65*se_women_`y_var'
}

replace p = _n
gen label = ""
replace label = "Original" if pos == 1 
replace label = "Prov-Year FE controlled" if pos == 2
replace label = "# of siblings controlled" if pos == 3 
replace label = "Sex ratio at birth controlled" if pos == 4
replace label = "GDP /capita at age 15 controlled" if pos == 5 
replace label = "Education resources controlled" if pos == 6 
replace label = "CSLs controlled"  if pos == 7
replace label = "All above" if pos == 8 

replace p = -p
foreach gender in "women" {
tw (scatter p coef_`gender'_n_birth , xscale(alt) xline(0) xtit("Coef. on Fines at 6-20 / Std", size(small))  mlabel(label) mlabp(12)  xlabel(-0.6(0.2)0.4, grid )) ///
(rcap up_`gender'_n_birth low_`gender'_n_birth p  , ytit(" ") horizontal lp(dash) lc(green)  legend(pos(5) ring(0)  size(small) lab(1 "Number of births") lab(2 "90% CI") col(1))  ylabel(-9 " " 0 " ") )
gr export "$path4/fig_c4e.eps",replace

}

use "$path2/child_gender_mort", clear 
gen wt = 1 
replace wt = 4 if year == 2005
keep if han_p == 1
replace die = die*100 
su die [aw = wt]
foreach y_var in "die" {
su `y_var' [aw = wt] if  han
sca std_child_`y_var' = `r(sd)'
}
sca list
use "$path3/coef_child_birth_control", clear
foreach y_var in "die"{
replace coef_child_`y_var' = coef_child_`y_var'/std_child_`y_var'
replace se_child_`y_var' = se_child_`y_var'/std_child_`y_var'
gen up_child_`y_var' = coef_child_`y_var' + 1.65*se_child_`y_var'
gen low_child_`y_var' = coef_child_`y_var' - 1.65*se_child_`y_var'
}
drop if p == 8 
replace p = _n
gen label = ""
replace label = "Original" if pos == 1 
replace label = "Prov-Year FE controlled" if pos == 2
replace label = "# of siblings controlled" if pos == 3 
replace label = "Sex ratio at birth controlled" if pos == 4
replace label = "GDP /capita at age 15 controlled" if pos == 5 
replace label = "Education resources controlled" if pos == 6 
replace label = "CSLs controlled"  if pos == 7
replace label = "All above" if pos == 8 

replace p = -p
foreach gender in "child" {
tw (scatter p coef_`gender'_die , xscale(alt) xline(0) xtit("Coef. on Fines at 6-20 / Std", size(small))  mlabel(label) mlabp(12)  xlabel(-0.1(0.02)0.04, grid )) ///
(rcap up_`gender'_die  low_`gender'_die  p  , ytit(" ") horizontal lp(dash) lc(green)  legend(pos(5) ring(0)  size(small) lab(1 "Child mortality") lab(2 "90% CI") col(1))  ylabel(-9 " " 0 " ") )
gr export "$path4/fig_c4f.eps",replace

}
 
 log close 

